home *** CD-ROM | disk | FTP | other *** search
/ Ian & Stuart's Australian Mac 1993 September / September 93.iso / Archives / Applications / Calculators / NumberCrunch 1.41 / Number Crunch II Demo / Library Routines / Text Examples / Random Numbers and Noise < prev    next >
Encoding:
Text File  |  1992-09-03  |  7.3 KB  |  258 lines  |  [TEXT/NCII]
  1. ##############
  2. # Random Numbers and Noise
  3. # 
  4. #  program InitRandom(seed)  # Intitializes all random routines; seed ~ 1.0 . 
  5. #   function Random # Returns a uniform random number between 0 and 1.
  6. #   function RandomArray(N) # Returns an array of N randoms.
  7. #   function RandomInRange(J,L) # Returns a random integer K, J<=K<=L
  8. # 
  9. #   function WhiteNoise(n) # Returns n-array of uniform, sigma=1 white noise.
  10. #   function GaussianNoise(N) # Returns array of gaussian noise, sigma=1.
  11. #  program MakeDistribution(Y,max,N,  x,P) # Bin randoms y into a distribution P
  12. # 
  13. #  This text file explains and gives examples of
  14. #  some of the routines in the file All Library Routines, which
  15. #  should be Opened before trying any of these examples.
  16. ##################
  17.  
  18. ##############
  19. # The xCOD itself is
  20. xRANDOM
  21.   # xCO2 xRANDOM(zInit,N,x:extended; RandStore:RealArray) , sets  0<x[1…N]<1 to random numbers. RandStore[1..100] must be preserved between calls. zInit <0 is a signal to initialize RandStore with zInit seed ; an external program.
  22.  
  23. ############
  24. # The interface routines.
  25. # Examples are given after each.
  26.  
  27. # This must be called once, initially, with a 0-10 or so
  28. # random seed.  The storage array RandStore must
  29. # exist between successive calls to xRANDOM.
  30.  program InitRandom(seed)  # Intitializes all random routines; seed ~ 1.0 . 
  31. .   var x
  32. .   #  Input:  seed = real random number ~ 1.0
  33. .   #  Output:  none.  Sets up RandStore.
  34. .   begin
  35. .     RandStore[1…100] = 0
  36. .     x[1…2] = 0
  37. .     xRANDOM(-314159*seed, size(x),x,RandStore)
  38. .   end
  39. .   
  40.  
  41. # If you want a different set of 'randoms', use a different
  42. # seed than the one here.
  43. InitRandom(0.133241)
  44.  
  45. ###########################
  46. # Here a number of different uses of xRandom.
  47. # Examples are given for each.
  48.  
  49. # This returns a single random number
  50. #  uniformly distributed from 0 to 1.
  51. #
  52.   function Random # Returns a uniform random number between 0 and 1.
  53. .   var x
  54. . #  Input: none
  55. . #  Ouput:   0 < Random < 1
  56. .   begin
  57. .     x=0
  58. .     xRandom(1,1,x,RandStore)
  59. .     Random=x
  60. .   end
  61. .   
  62.  
  63. # For example, 
  64. Random               # type this, hit return
  65.   0.70898038     # and get this the 1st time
  66. Random
  67.   0.57992348     # but this the second.
  68.  
  69.  
  70. # Returns an array of N random numbers between 0 and 1.
  71. # Much faster than successive calls to Random, above.
  72. #
  73.   function RandomArray(N) # Returns an array of N randoms.
  74. .   var x
  75. .  #   Input:    N = positive integer
  76. .  #   Output:  RandomArray = x[1…N] random numbers, all 0 -> 1.
  77. .   begin
  78. .     x[1…n]=0
  79. .     xRandom(1,n,x,RandStore)
  80. .     RandomArray = x
  81. .   end
  82. .   
  83.  
  84. # For example,
  85. RandomArray(5)
  86.   [0.85127, 0.2307, 0.65812, 0.64353, 0.11952] 
  87.  
  88.  
  89. # Gives a random integer K in the range  J <= K <= L.  
  90. #  Note that the Int(x) function rounds to
  91. #  the nearest even integer.
  92. #
  93.   function RandomInRange(J,L) # Returns a random integer K, J<=K<=L
  94. .   var x, a, b
  95. . #   Input:   J, L = integers
  96. . #   Ouput:   RandomInRange = random integer from J to L.
  97. .   begin
  98. .     a = J
  99. .     b = L
  100. .     a=a-.5+1e-17
  101. .     b=b+.5-1e-17
  102. .     x=0
  103. .     xRandom(1,1,x,RandStore)
  104. .     x=a+x*(b-a)
  105. .     RandomInRange=Int(x)
  106. .   end
  107. .   
  108.  
  109.  
  110. # For example, to simulate the roll of a die,
  111. RandomInRange(1,6)
  112.   4 
  113. RandomInRange(1,6)
  114.   4 
  115. RandomInRange(1,6)
  116.   2 
  117.  
  118. # Here are several more trials:
  119. for j=1:20 do d[j]=RandomInRange(1,6)
  120.   d[1…20] = [3, 5, 1, 2, 5, 3, 5, 4, 5, 6, 4, 4, 6, 6, 2, 3, 4, 2, 6, 5] 
  121.  
  122.  
  123. # Here are two that calculate some statistics.
  124. # These are also in the file Miscellaneous Routines.
  125.   Mean(V) = sum(V)/size(V) # Returns mean of an array V.;
  126.   StandardDeviation(V) = sqrt( Mean(V^2) - Mean(V)^2);
  127.  
  128. # For example, 
  129. Mean(d)   # where d is the dice rolls done above
  130.   4.05 
  131. StandardDeviation(d)
  132.   1.49916644 
  133.  
  134. # The average of a number of dice rolls should be (1+6)/2 = 3.5.
  135. # In N trials, the error in the average should be about 1/sqrt(N).
  136. # For d, with 20 trials, the expected error is then
  137. 1/sqrt(20)
  138.   0.2236068 
  139. # while the actual error actually was
  140. 4.05 - 3.5
  141.   0.55 
  142. # which means that usually we'll get a bit closer than 4.05.
  143.  
  144. # Another trial gave me
  145. for j=1:20 do d[j]=RandomInRange(1,6)
  146.   d[1…20] = [2, 6, 6, 5, 4, 1, 6, 5, 1, 5, 4, 3, 1, 1, 5, 4, 3, 4, 6, 5] 
  147. Average(d)
  148.   3.85 
  149. # which is closer.
  150.  
  151.  
  152. # The average of Random itself is clearly 0.5, halfway from 0 to 1.  Since the average of x^2 is integral(x^2 from 0 to1) = 1/3, 
  153. # the StandardDeviation is
  154. SigmaRand = sqrt(1/3-1/4)
  155.   SigmaRand = 0.28867513 
  156. AverageRand = 0.5
  157.  
  158. # With this we can define a white-noise routine, which
  159. # generates an N component array of  0 mean, 
  160. # standard dev.=1 numbers.
  161. #
  162.   function WhiteNoise(n) # Returns n-array of uniform, sigma=1 white noise.
  163. .   var AverageRand, SigmaRand
  164. .   # Input: N = positive integer
  165. .   # Output: WhiteNoise = array W[1…N] of 0 mean, variance = 1,
  166. .   #                                 uniform random numbers.
  167. .   begin
  168. .     AverageRand = 0.5;  SigmaRand = sqrt(1/3-1/4);
  169. .     WhiteNoise = (RandomArray(n)-AverageRand)/SigmaRand;
  170. .   end
  171. .   
  172.  
  173. # Some white noise.
  174. w = whiteNoise(1000);
  175.  
  176. # And its statistics.
  177. Mean(w)
  178.   -0.04581557           # This should be about zero.
  179. StandardDeviation(w)
  180.   0.97346018             # This should be about 1.
  181. 1/sqrt(1000)
  182.   0.03162278             # within about this tolerance.
  183.  
  184. # Plot [1…100, WhiteNoise(100)] to see a graph of this.
  185. # Each time you hit the plot button it should change.
  186.  
  187.  
  188. # Finally, noise with a gaussian spectrum. 
  189. # For an explanation of this, see a text like
  190. # Numerical Recipes by W. Press et al.
  191. # This returns an N component array of 0 mean
  192. # standard deviation=1 Gaussian (bell-curve)
  193. # distributed noise. 
  194. #
  195.   function GaussianNoise(N) # Returns array of gaussian noise, sigma=1.
  196. .   var x,y
  197. . #  Input:   N = positive integer
  198. . #  Output:   GaussianNoise = G[1…n] = array of 0 mean, sigma=1 randoms.
  199. .   begin
  200. .     x = RandomArray(n)
  201. .     y = RandomArray(n)
  202. .     GaussianNoise = sqrt(-2 ln(x)) sin(2πy)
  203. .   end
  204. .   
  205.  
  206. g = GaussianNoise(300);
  207.  
  208. # To see a scatter plot of this, just plot yArray=g in Quick Setup 
  209. # in the graph.
  210.  
  211. # Seeing the actual probability distrubution is trickier, since you
  212. # want to essentially how many points of g lie in each dx range.
  213. # This counting can be done by using the g itself as an index into
  214. # another array.
  215.  
  216. # This calculates the probability distribution P of Y,
  217. # calculated by binning it into N+1 boxes
  218. # from -max to max. x is also defined, so that plotting [x,P] 
  219. # shows the distribution
  220. #
  221.  program MakeDistribution(Y,max,N,  x,P) # Bin randoms y into a distribution P.
  222. .   var min,max, yScal,j, Ny
  223. .  #  Input:  Y[1…j] = array of random numbers
  224. .  #             max = range of Y to examine, from y=-max to max
  225. .  #             N = number of bins to count Y in.
  226. .  #  Output:
  227. .  #             P[1…(n+1)] distribution
  228. .  #             x[1…(n+1)] points where P is defined.
  229. .   begin
  230. .     P[-N/2…N/2] = 0;   
  231. .     yScal=y*(N/2)/max
  232. .     x=-N/2…N/2;  x =x*(max/(N/2))
  233. .     Ny=size(Y)
  234. .     j=1
  235. .     while j<Ny do
  236. .       begin
  237. .         P[yScal[j]] = P[yScal[j]]+1
  238. .         j=j+1
  239. .       end
  240. .     P = P/sum(P)
  241. .   end
  242. .   
  243.  
  244.  
  245. # In the case of the gausian "g" above, it ranges from
  246. minimum(g)
  247.   -2.58766783 
  248. Maximum(g)
  249.   3.3347715 
  250. # so if we bin it in, say, 30 bins (want some of 100 in each),
  251. # then [x,P] are set by:  (This is slow, give it a few minutes.)
  252. MakeDistribution(g,3.4,30, x,P)
  253.  
  254. # Plotting [x,P] now shows (roughly) a gaussian bell curve. Increase
  255. # the array sizes for a nicer picture.
  256.  
  257.  
  258.